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We have numerically studied the instability of the spherically symmetric 
Qh| standing accretion shock wave against non-spherical perturbations. We have in 

O ' mind the application to the collapse-driven supernovae in the post bounce phase, 

^ \ where the prompt shock wave generated by core bounce is commonly stalled. 

We take an experimental stand point in this paper. Using spherically symmetric, 
completely steady, shocked accretion flows as unperturbed states, we have clearly 
^ ' observed both the linear growth and the subsequent nonlinear saturation of the 

d \ instability. In so doing, we have employed a reahstic equation of state together 

with heating and cooling via neutrino reactions with nucleons. We have done a 
mode analysis based on the spherical harmonics decomposition and found that 
the modes with £ = 1,2 are dominant not only in the linear regime, but also after 
the nonlinear couplings generate various modes and the saturation occurs. Vary- 
ing the neutrino luminosity, we have constructed the unperturbed states both 
with and without a negative entropy- gradient. We have found that in both cases 
the growth of the instability is similar, suggesting the convection does not play a 
dominant role, which also appears to be supported by the recent linear analysis 
of the convection in accretion flows by Foglizzo et al. The real part of the eigen 
frequency seems to be mainly determined by the advection time rather than by 
the sound-crossing time. Whatever the cause may be, the instability is favorable 
for the shock revival. 
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Subject headings: supernovae: collapse — neutrinos — hydrodynamics — insta- 
bility 

1. Introduction 

Over the last decades, it has been observationally suggested that the core-collapse su- 
pernovae are generally aspherical (Wang et al. 1996, 2001, 2002). The most unequivocal 
example is SN1987A. Recent HST images of SN1987A are directly showing that the expand- 
ing envelope is indeed elliptical with the long axis aligned with the rotation axis inferred 
from the ring. The aspect ratio and the position angle of the symmetry axis are consistent 
with those predicted earlier from the observations of speckle and linear polarization. What is 
more, the linear polarization became greater as the time passed (Wang et al. 2001; Leonard 
et al. 2001), the fact which has been used to argue that the central engine of the explosion is 
responsible for the non-sphericity (Khokhlov et al. 1999; Wheeler et al. 2000). Rather com- 
mon detections of the linear polarization in other core-collapse supernovae seem to suggest, 
as mentioned at the beginning, that the core-collapse supernovae are globally asymmetric in 
general. 

Various physical ingredients have been proposed as a possible cause of the asymmetry 
so far: convections below the shock wave as well as inside the proto neutron star (Herant 
et al. (1994); Burrows et al. (1995); Janka & Mueller (1996), see for a review Janka et al. 
(2001)), density inhomogeneities produced prior to core-collapse (Burrows & Hayes 1996; 
Fryer 2004), rapid rotation of the core (Monchmeyer & Miiller 1989; Yamada & Sato 1994; 
Shimizu et al. 2001; Kotake et al. 2003; Yamasaki & Yamada 2005; Walder et al. 2005), 
and magnetic fields (Ardeljan et al. 2000; Akiyama et al. 2003; Kotake et al. 2004; Yamada 
& Sawai 2004; Kotake et al. 2005) (see for a review, Kotake et al. (2005) and references 
therein). These studies have also been motivated by the interest in the consequences of 
the asymmetric motions to the explosion mechanism itself, since it has been expected that 
the deviation from the spherical symmetry will be helpful for the shock revival one way or 
another. 

Recently, Blondin et al. (2003) demonstrated numerically that yet another hydrody- 
namical mechanism may be conspiring to drive non-spherical motions in the flow below the 
shock wave. The so-called standing accretion shock instability (SASI) is supposed to be a 
non-local hydrodynamical instability possibly caused by the cycle of the inward advection of 
velocity- and entropy-perturbations and outward propagation of acoustic waves, with fluc- 
tuations amplified after each cycle. This mechanism of SASI was originally studied in linear 
analysis by Foglizzo (2001, 2002) in the context of accreting black holes (see also Houck & 
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Chevalier (1992)). Adding small non-spherical perturbations to the spherically symmetric, 
isentropic, steady, post-shock accretion flows, Blondin et al. (2003) found in their numerical 
simulations that the perturbations grow up to the nonlinear regime with clear dominance of 
£ = 1 at first and 1 = 2 later, leading to the global deformation of the shock wave. Here I 
stands for the azimuthal index of the Legendre polynomials. One lesson to learn is that we 
should not impose the symmetry with respect to the equatorial plane in the simulations. 

As mentioned already, since the large deviation from the spherical symmetry may have 
an important consequence to the explosion itself, the finding of Blondin et al. (2003) has 
much interest of other researchers. In their first paper (Blondin et al. 2003), the neutrino 
heating and cooling are entirely ignored and the flow is assumed to be isentropic. In the 
recent paper (Blondin & Mezzacappa 2005), the authors took into account the cooling term 
of a simple analytic form just as in Houck & Chevaher (1992), but no heating included 
yet. On the other hand, Scheck et al. (2004) demonstrated that similar asymmetric motions 
with no equatorial symmetry occur in their most realistic mimcrical models. Although their 
results show that the neutrino processes will not nullify SASI, the growths and saturations of 
individual modes under neutrino-irradiation are not clear, since they used highly complicated 
flows as an underlying model. 

Our standing point in this paper is somewhere in between these works. We study SASI 
by 2D axisymmetric hydrodynamical simulations. Although we have in mind the application 
to the supernova core in the shock-stagnation phase, we take an experimental stance as in 
Blondin et al. (2003). On the one hand, we employ a realistic equation of state (Shcn et 
al. 1998) and take into account the heating and cooling of matter via neutrino emissions 
and absorptions on nucleons. As an underlying model, on the other hand, we utilize the 
spherically symmetric, steady, shocked accretion flows (Yamasaki & Yamada 2005), which 
is stable against radial perturbations. Although this is certainly a crude approximation to 
what we found in the realistic simulations, it will enable us to do clear mode-analyses from 
the linear growths through the nonlinear couplings among various modes up to the eventual 
saturation of SASI. Due to the neutrino-heating, some initial models have a convectively 
unstable region in the classical sense (see Foglizzo et al. (2005)) and SASI is inevitably 
mixed with convection. By lowering the neutrino luminosity, however, we can also construct 
models with no convectively unstable region. By comparing these models, we can assess the 
relative strength of these instabihties. We also discuss the implications that SASI might 
have for the shock revival. 

The plan of this paper is as follows. We describe the numerical methods in section 2. 
The main numerical results are shown in section 3. We conclude this paper in section 4. 
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2. Numerical Methods and Models 

2.1. Basic Equations 

Assuming the axisymmetry of the system in this paper, we numerically study the growth 
of the non-spherical instability in the accretion flow through the shock wave onto the pro- 
toneutron star. The unperturbed steady accretion flows and the shock waves are assumed 
to be spherically symmetric. We take into account the heating and cooling of accreting 
matter via neutrino absorptions and emissions by free nucleons. Only the region outside the 
neutrino sphere is considered. 

The basic evolution equations are written as follows, 

^ + pV-v^O, (1) 

= -VP - pV*, (2) 



dt 



Qn, (4) 



$ = -^, (5) 

r 

where p, v, e, P, Y^, $ are density, velocity, internal energy, pressure, electron fraction, and 
gravitational potential of the central object, respectively. The self-gravity of matter in the 
accretion flow is ignored. Qe and Qn are related with the interactions with neutrinos and 
are explained in more detail in the next section. We denote the Lagrangian derivative as 
d/dt and r is the radius. 

The numerical code for hydrodynamic computations employed in this paper is based 
on the ZEUS-2D (Stone & Norman 1992), which is an Eulerian code based on the finite- 
difference method and employs an artificial viscosity of von Neumann and Richtmyer type 
to capture shocks. We have made several major changes to the base code to include the 
microphysics as described in the following sections. First, we have added the equation for 
electron fraction (Eq. (4)), which is solved in the operator-splitting fashion. Second, we have 
incorporated the tabulated realistic equation of state (EOS) based on the relativistic mean 
field theory (Shen et al. 1998) instead of the ideal gas EOS assumed in the original code. 

Spherical coordinates are used. No cqTiatorial symmetry is assumed and the computation 
domain covers the whole meridian section with 60 angular mesh points, except for a model 
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in which we adopted 120 angular grid points. Since the latter model did not produce any 
significant difference from other models, we will report in the following the results obtained 
from the models with 60 angular mesh points. We use 300 radial mesh points to cover 
^ r < Tout = 2000 km, where Tin is the inner boundary and chosen to be the radius of 
neutrino sphere, r^, defined later. 



2.2. Neutrino Absorptions and Emissions by Free Nucleons 

Qe and Qn in Eqs. (3) and (4) are concerning the absorptions and emissions of electron- 
type neutrino and anti-neutrino z/g on free nucleons. We summarize here the expressions 
for these terms together with the approximations we employ in this paper. 

The neutrino absorptivity on free neutron (i/g+n — > e~-|-p) is written as follows (Bruenn 
1985; Rampp & Janka 2002), 



e) = Qn.vi'igl + 9^) [1 - Fe(e + A)] (e + A) V(e + A)^ - m^c^, (6) 



where me, c, e are electron mass, speed of light, and neutrino energy, respectively. Q is 
defined as follows. 
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where the characteristic cross section of weak interaction, ctq, is given as 

. ^(^;f^/)^ ^ 1.761 . 10-- cm^ (8) 

In the above equation, Gp is the Fermi constant. A is the difference of rest mass energies of 
neutron and proton; 

A = TTijiC^ — mpC^. (9) 

The nuclcon form factors for the vector and axial vector currents are set to be ^fy = 1 and 
gA = 1.23, respectively. The degeneracy factor, ?7np, is defined by 

Fj is the Fermi-Dirac distribution function for species i; 

F,(e,) = 1/ [1 + exp((e, - iii)/kBT)] , (11) 

where /cb, T, /li are Boltzmann constant, temperature, and chemical potential including rest 
mass, respectively. The emissivity of can be obtained from the absorptivity by detailed 
balance as 

j{e) = exp [- (e - (//p - Hn + He)) /kBT] ^(e). (12) 
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The absorptivity of anti- neutrino on free proton (i^e+P e++n) is similarly represented 

by 

= SVpni^gl + 9y) [1 - ^e+(e - A)] (e - A) V(e - A)^ - ^^(^6(6 - A - m,c'), (13) 
where 0(e) is the step function; 

The degeneracy factor, r/pn, is defined by exchanging the contributions from proton and 
neutron in Eq. (10). The emissivity is obtained again by the detailed balance just like for 
Eq. (12) together with the relation of (1^+ = —/le, 

j{e) = exp [- (e - (//n - A*p - A*e)) /kBT] K{e). (15) 

Using the absorptivities and emissivities given above, Qe and Qn can be expressed as 
follows, 

Qe = -J^s I ^'de m - (j(e) + n{e)) f{r, e)] , (16) 

= V (sS^ r ''"^^ ^^'^ " ^^^'^ ^ "^'^^ ' ^^^^ 

i — —1 (for electron-type neutrino) 
i — 1 (for electron-type antineutrino) ' 

where /(r, e) is the distribution function of neutrino. These rates are calculated for i/^ and 
Pe separately and summed for the source terms in Eqs. (3) and (4). 

Since we deal with the optically thin region outside the neutrino sphere in this paper, 
we do not have to solve the transport equation for neutrinos. We assume that the neutrino 
distribution functions are approximated by the Fermi-Dirac distribution with a vanishing 
chemical potential; 

where the geometrical factor is taken into account for normalization. Note that although 
the angular dependence in the above distribution function is assumed to be isotropic, it is 
entirely irrelevant for the absorption and emission rates. We further assume for simplicity 
in the following that the luminosity Lj,, temperature Ti, and neutrino sphere rj, are related 
by the following equation, 

= ^aT^ ■ Anrl (19) 
where a is the Stefan-Boltzmann constant. 
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2.3. Initial and Boundary Conditions 

In order to clearly see the linear growths of the instability, it is critically important to use 
a well-defined unperturbed state as an initial condition for simulations. For this purpose, 
we employ the spherically symmetric steady accretion fiows through the standing shock 
wave in this paper. Following Yamasaki & Yamada (2005), we solve the time-independent 
hydrodynamical equations from the shock front down to the inner boundary, 



Anr'^pvr = M, (20) 

dvr 1 dP GM„ 

dr p dr 



-r^ + -- + ^-^^ (21) 



de P dp Qe 

Vr^ = , (22) 

dr p^ dr p 
dY^ 

^''~dr ^ ^'^^^ 
where M and Vr are the mass accretion rate and radial velocity, respectively. Given the 
post-shock values for density, radial velocity, entropy and electron fraction (ppost) f^post) •^post; 
y^post), the shock radius is determined so that the density obtained at the inner boundary 
should agree with the fixed value, pin = 10^^ g/cm^. 

The post-shock values are calculated from the corresponding pre-shock values by the 
Rankine-Hugoniot relations. The pre-shock values, on the other hand, are obtained from the 
mass accretion rate, the outer boundary conditions for entropy and electron fraction (sout) 
y^out) the assumption that matter falls freely with no interaction with neutrinos outside 
the shock; 

Ve = -yj2GMjR,, (24) 

*pre — *out) (26) 
■^pre ~ ^out) (27) 

where Rg is the shock radius. The entropy and the electron fraction at the outer boundary 
are assumed to be Sout — 3 /ce/baryon and ygout — 0-5, respectively. 

In the numerical simulations, axisymmetry but no equatorial symmetry is assumed. At 
the outer boundary, we adopt the fixed boundary condition consistent with the unperturbed 
state. On the other hand, the free-fiow-in boundary condition is used at the inner boundary. 
A reahstic tabulated EOS by Shen et al. (1998) is used both for the simulations and for the 
preparation of the initial conditions. 
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For all the models investigated in this paper, the mass accretion rate and the mass of 
central object are fixed to be M = 1 Mq/s and Afin = 1.4 Mq, respectively. The neutrino 
temperatures are also fixed to T^^ = 4 MeV and Tp^ = 5 MeV, the typical values in the post- 
bounce phase. On the other hand, we systematically vary the neutrino luminosity within the 
range of L,y = 3.0-7.0 ■ 10^^ erg/s for different models. Note that these values are constant 
in time for each model. This is necessary to realize the steady unperturbed states. 

To study the instability, we add angular-dependent perturbations, Svr{r, 9), to the radial 
velocity initially; 

y^(^r,e) =vl^{r) + 6vr{r,e), (28) 

where vj.^ is the spherically symmetric unperturbed velocity. The i = 1 {Svr{r,6) oc cos^) 
single-mode perturbation and the random multi-mode perturbation are investigated. In the 
former case, the perturbation amplitude is set to be 1% of the unperturbed velocity. For 
the latter case, the amplitude is less than 1% for each radial direction. In both cases, the 
assumed initial amplitudes are found to be small enough to observe clearly the linear growths 
of the instabihty. 



3. Results 

3.1. Basic Features 

The unperturbed spherically symmetric steady accretion flows are displayed in Fig. 1 
for Lj, — 3.0,5.5,6.0 • 10^^ erg/s. As L^, increases, the width between the shock and the 
neutrino sphere, Ws = Rs — fi,, becomes larger with both the neutrino sphere and the shock 
radius getting larger. In the cases of L,y = 5.5,6.0 ■ 10^^ erg/s, the region with a negative 
entropy-gradient is formed behind the shock, where the net heating rate is positive. The 
supernova core is convectively unstable in this region. In the case of L^, = 3.0 ■ 10"''^ erg/s, on 
the other hand, the net heating rate is always negative over the whole region from the shock 
down to the neutrino sphere, since the shock radius is smaller and the temperature is higher, 
implying more efficient coolings. As a result, there is no region with a negative entropy- 
gradient. Hence we can see in this case the instability solely caused by SASI. It is also noted 
that this model should be similar to the models considered by Blondin & Mezzacappa (2005) 
with only simple cooling terms taken into account. 

We use the profiles in Fig. 1 as initial conditions for hydrodynamical simulations in 
2D. Since they are obtained by solving Eqs. (20)-(23), a shght inconsistency is inevitably 
introduced by the mapping. For example, the shock is smoothed over a few computational 
grid points with the artificial viscosity in the simulations while the initial conditions satisfy 
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the Rankine-Hugoniot relation across the infinitesimal distance. As a result, transient oscil- 
lations are commonly found before they are damped away. In fact, we have done spherically 
symmetric computations for the unperturbed initial conditions using the same dynamical 
code before doing 2D simulations with perturbations. The purpose is two-fold. First, we can 
obtain this way "the steady states" for the dynamical code, which, as we have emphasized 
repeatedly, is critically important to see the linear growth of the instability. These simu- 
lations are also meant for the confirmation of the stability of the initial conditions against 
radial perturbations predicted by Yamasaki & Yamada (2005). 

The temporal evolutions of the shock radii in the spherically symmetric simulations are 
shown in Fig. 2 for the three cases in Fig. 1 together with for L^j = 6.5,7.0 ■ 10^^ erg/s. 
The oscillations of shock radii are particularly clear in the early phase. The oscillation 
frequency is found to be inversely proportional to the neutrino luminosity. For the neutrino 
luminosity less than 7.0 • 10^^ erg/s, the oscillations are gradually damped and the shock 
approaches to the equilibrium radius, which is a bit different from the initial value. In 
the case of L^, = 7.0 • 10^^ erg/s, however, the oscillation is slowly amplified and appears 
to be never settled to equilibrium. It is noted that this model is expected to be stable 
according to Yamasaki & Yamada (2005). The apparent discrepancy may be ascribed to 
the numerical errors inherent to the dynamical simulations as mentioned above. It should 
be also mentioned, however, that the large-amplitude oscillations as observed in the figure 
may cause instability even for linearly stable configurations. Since the main purpose of this 
paper is to investigate the stability for non-radial perturbations, we focus in the following 
only on the models with L^, < 6.0 • 10^^ erg/s, in which the radial oscillations are damped 
within ~ 100 ms. After these transient feateures are sufficiently settled, we add non-radial 
velocity perturbations explained in the previous section to the profiles obtained this way and 
start 2D simulations. 

Figure 3 shows in the meridian section the distributions of entropy (the left half of 
the panel) and density (the right half) for the models with = 5.5,6.0 ■ 10^^ erg/s after 
1% of the £ = 1 single-mode velocity perturbation is added. For both models, we observe 
the growth of the perturbations. In the case of L^, = 5.5 ■ 10^^ erg/s, the shock surface is 
deformed at first by the increasing amplitude of the non-radial mode and then begins to 
oscillate with a large amplitude. In the case of — 6.0 • 10^^ erg/s (right panels), on the 
other hand, in addition to the oscillations of the shock surface, we observe the substantial 
increase of the average shock radius as the time passes. In fact, after t — 400 ms, the shock 
radius continues to increase and appears to produce an explosion. Since the model is stable 
against radial perturbations as mentioned above, the non-radial instability and the neutrino 
heating therein are responsible for the explosion. We think that this is a reconfirmation of 
the claim that the instability, whatever the cause, behind the shock is helpful for the shock 
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revival. 



In order to analyze further in detail the evolution of the instability, we conduct a mode 
analysis as follows. The deformation of the shock surface is decomposed into the spherical 
harmonic components; 



Since the system is axisymmetric, only m — harmonics, nothing but Legendre polynomi- 
als, show up. The coefficients, a^, can be calculated by the orthogonality of the Legendre 
polynomials; 



The position of the shock surface, Rs{0), is estimated from the iso-entropic surface of s = 5 
in the following analyses. 

The temporal evolutions of the average shock radius are obtained from the i = com- 
ponent in the above decomposition and \ ARs/ Rsfl\ = \ {ao ~ -Rs,o)/-Rs,o| ^''^^ plotted in Fig. 4, 
where i?s,o is the initial shock radius. As can be seen from the figure, after the shock sur- 
face oscillates for about 100 ms as in the symmetric simulations, it begins to increase. For 
the models with L^, = 3.0,5.5 • 10^^ erg/s, the shock surface is settled in ~ 50 ms to a 
quasi-steady state with a radius ~ 10% larger than the initial value, and keeps oscillating 
around it thereafter. In the case of L^, = 6.0 • 10^^ erg/s, on the other hand, a continuous 
increase of the average shock radius is found, which will eventually lead to the explosion men- 
tioned above. It is emphasized that even for the model with no negative entropy-gradient 
{Li, = 3.0 ■ 10^^ erg/s), the exponential increase of the shock radius occurs during the same 
period. This demonstrates clearly that the instability has nothing to with the convection at 
least for this model. The relative contributions of SASI and convection for other models will 
be discussed in the next section. 



In this section, we discuss both the linear and nonlinear growths of the instability 
in greater detail. As reported in Foglizzo (2001, 2002); Blondin ct al. (2003), the most 
remarkable feature of SASI is the dominance of the modes with ^ = 1,2 in the instability. 
Hence we start with the analysis of the models with £ = 1 single-mode velocity perturbation 
imposed initially. 




(29) 




(30) 



3.2. Linear and Nonlinear Growths of Instability 



Fig. 5 shows the temporal evolutions of the £ = 1 and £ = 2 modes. Although the result 
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for L,^ — 5.5 • 10^^ erg/s is plotted as a reference case, the qualitative feature is common 
to other models. As can be clearly seen, the evolution is divided into two phases. The 
initial phase lasting for ~ 100 ms represents the linear phase, where the amplitude of each 
mode grows exponentially. It is noted that even for the single-mode perturbation, the second 
harmonics is generated by the nonlinear coupling and grows exponentially also. In the linear 
phase, however, the second harmonics is always much smaller than the fundamental i = 1 
mode. After ~ 100 ms, the amplitude of the fundamental mode is saturated at ~ 10% level 
and the nonlinear phase starts. It is found that soon after the nonlinear phase begins, the 
amplitude of the second harmonics becomes comparable to the fundamental mode. This 
corresponds to the instability found in the simulations by Blondin et al. (2003). 

In order to obtain the linear growth rates for the fundamental mode = 1), we fit the 
amplitude of the mode, ai{t), with the following expression. 



where 7 and 00 are the exponential growth rate and the characteristic oscillation frequency 
in the linear phase, respectively. The least square fitting is done for these parameters and 
the overall normalization. A, and the phase, 5. The dash dotted line in Fig. 5 represents the 
result. The obtained values of 7 and uj are listed in Table 1 for the models with different L^,. 

According to Foglizzo (2001, 2002), the instability is produced by the cycle of the in- 
ward advection of the velocity and entropy fiuctuations and the outward propagation of the 
pressure fiuctuations (sec also Blondin et al. (2003)). If this is correct, the characteristic 
oscillation frequency refiects the the cycle period. 



where Cg is the sound velocity. The values of a; estimated this way are also given in Table 1. 
They are found to agree quite well with the numerical results. This is in contrast with the 
recent claim by Blondin & Mezzacappa (2005) that the propagation of pressure perturbations 
is responsible for the instability. Their simplified treatment of the cooling terms may be 
responsible for the difference. Our simulations indicate clearly that the time scale associated 
with the advection plays an important role in the instability. We can obviously find the 
Lj^-dependence of the oscillation frequency, cu, which can be understood as follows. As also 
presented in Table 1, the width between the shock and the neutrino sphere, Wg — Rs,equii — T^i^, 
becomes larger as the luminosity increases. As a result, the oscillation frequency is expected 
to become lower since the fiuctuations traverse longer distances in the cycle. 



ai (t) = A exp(7t) sin(u;t + S) 



(31) 




(32) 



The growth rate, 7, on the other hand, has little dependence on the luminosity as 
shown in Table 1, As mentioned already, the instability occurs even in the case of L^, — 
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3.0 • 10^^ erg/s, where the negative entropy-gradient is not formed and the structure is stable 
for convection. The interesting thing here is the fact that the hnear growth rates are not so 
different between the models with and without a negative entropy-gradient. This suggests 
that SASl plays a dominant role in driving the non-radial motions even when the convection 
is also expected to occur. This issue will be further discussed later in this section again. 

The simulations discussed so far have been done with the initial velocity perturbation 
including only the fundamental mode {i = 1). However, other modes with i >2 also develop 
rapidly from ~ 50 ms. After ~ 100 ms, the amplitude of £ = 2 mode becomes the same 
order as that of the fundamental mode, which has already been saturated by this time. This 
marks the beginning of the nonlinear phase. In fact, the £ = 2 mode is also soon saturated. 
This transition from linear to nonlinear phase corresponds to the time of the rapid increase 
of the average shock radius shown in Fig. 4. Since the average radius is nothing but the 
£ — mode, this can be interpreted as a result of the nonlinear coupling of this mode with 
the fundamental mode and the ensuing saturation. As will be mentioned later, since the 
expansion of the shock is crucial to trigger the shock revival and eventual explosion, it is 
important that the neutrino luminosity is sustained for ~ 100 ms, the typical saturation 
time for SASL 

Next we discuss the models with the random multi-mode velocity perturbations. So 
far we confirmed that the £ = 1 mode is indeed unstable to SASl. The models are meant 
to see which mode is dominant. In so doing, we also study the influence of the existence 
of a negative entropy- gradient. In Fig. 6, the temporal evolution of the spectrum of the 
spherical harmonics. The spherically symmetric component, the £ = mode, is omitted 
in the figure. The cases without and with a negative entropy-gradient are shown in the 
left and right panels, respectively. It is obvious that the modes with small £'s, especially 
those with i — 1,2, grow rapidly in the hnear regime {t < 100 ms). This is particularly 
the case for the model without a negative entropy-gradient {L^, — 3.0 • 10^^ erg/s) and the 
growths of the modes with £ > 10 are negligibly small. With a negative entropy-gradient 
{L,y = 5.5 ■ 10^^ erg/s), the broadening of spectra to larger £ modes is observed although the 
dominance of smaller £ modes can be still found. The convective instability may enhance 
the growth of higher harmonics in the linear phase. The similarity of the two cases suggest 
again that SASI is dominant over the convection even when the latter is operating. 

Recently Foglizzo et al. (2005) discussed the linear stability for convection in the accre- 
tion flows in the supernova core. They found that the classical criterion for convection, that 
is, the negative entropy-gradient is not sufficient for the accretion flows, since the limited time 
is available for growth. Although the classical convection has greater linear-growth rates for 
modes with larger wave numbers, they claimed that there are minimum (/cmin) and maximum 
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(^max) wave numbers for unstable modes in the accretion flows, and that the growth rates 
are also modified. The important parameter is the ratio of the advection time through the 
gain region divided by the local timescale of buoyancy, x, given in Eq. (40) in their paper. 
Applying their formula to our models, wc obtain % = 4 ... 7 for Lj^ = 5.5 ... 6.5 ■ 10^^ erg/s, 
with larger values for greater luminosities. Hence, the initial configurations are unstable 
against convection for these models, since the criterion, % > 3 is satisfied. The minimum 
and maximum wave numbers estimated in our models are /Cmin = 2...6 • 10~^ cm~^ and 
^max — 1 • 10~^ cm~^, respectively. The smaller k^i^^ corresponds to larger luminosities. 
They roughly correspond to the minimum and maximum indices in the spherical harmonics, 
^min = 2 ... 4 and ^max = 70 ... 85, respectively. The lower i^nm and larger ^^ax are obtained 
for higher luminosities. These number appear to be consistent with the spectrum shown 
in Fig. 6 as also inferred from Fig. 5 in Foglizzo et al. (2005). Since the classical growth 
rate of convection is comparable to that of SASl in our models and the true growth rate 
of convection will be much smaller than the classical estimation, we think that SASI is a 
dominant driving force for the non-radial motions we observed so far. 

It is also interesting to note that the modes with £=1,2 are dominant in the nonlinear 
regime, which begins after ~ 100 ms. As clearly seen in the broadening of the spectra in 
Fig. 6, various modes are amplified by nonlinear couplings with the dominant modes in this 
phase. The spectra are again broader for the model with a negative entropy- gradient. In 
both models, however, the dominance of the modes with £ = 1, 2 is remarkable. This should 
correspond to the large deformations of shock wave found in the numerical simulations by 
Blondin et al. (2003) and Scheck et al. (2004). In order to make clear the reason for the 
dominance of these modes in the nonlinear regime, it will be required to study the nonlinear 
couplings of various modes in more detail. 



3.3. Neutrino Heatings in SASI 

Finally, we discuss the influence of the instability on the heating of matter by neutrinos. 
Figure 7 shows the distributions of the net heating rates in the meridian section for the 
models with Li, = 5.5 • 10^^ erg/s (left panel) and L,, = 6.0 • 10^^ erg/s (right panel). In 
both cases, the shock wave is asymmetric with respect to the equator, the characteristics for 
SASI. For = 5.5 ■ 10^^ erg/s, however, the average shock radius is only slightly larger than 
that in the unperturbed state. As a result, both the heating and cooling region arc not much 
changed during the period from 100 ms to 200 ms. They are just oscillating. In the case of 
Lj, = 6.0 • 10^^ erg/s, on the other hand, not only the shock radius but also the heating region 
is getting larger as the time passes. It is further seen that the coohng region at 200 ms is 
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smaller than that at 100 ms. This is caused by the non-radial flows which carry down colder 
matter more efficiently than the spherical flows. The temperature near the inner boundary 
is lowered as a result. Since the cooling rates are roughly proportional to T^, the net cooling 
region is shrunk. This mechanism has been discussed in the context of convections over the 
years. It is applicable to SASl as it is. This broadening of the heating region and the shrink 
of the cooling region thrusts the shock wave outwards further, the positive feed back which 
eventually leads to the explosion as mentioned already. 

This positive feed back cycle is confirmed by the comparison of the temporal evolution 
of the shock radius displayed in Fig. 4 and that of the net heating rate integrated over the 
volume inside the shock wave shown in Fig. 8. We can see the gradual and continuous 
increase of the net heating after ~ 200 ms coincides with the expansion of the shock wave. 
This should be compared with the case for Li, — 5.5 • 10^^ erg/s, which does not produce any 
explosion. In this case, the average shock radius is settled to the value ~ 10% larger than 
the initial radius after ~ 150 ms while the net heating rate is fluctuated around a constant 
value during the same period. 

Figure 9 shows the temporal evolutions of various energies integrated over the region 
inside the shock wave for the model with L,^ — 6.0 • 10^^ erg/s. The explosion energy shown 
in the figure is defined as the total energy summed only over the fluid elements which have a 
positive total energy and a positive radial velocity simultaneously, as was done, for example, 
by Yamada & Sawai (2004). It is again clearly seen that the increase of the internal energy 
after ~ 200 ms is responsible for the shock revival. Note that the explosion energy remains 
zero for a while after the shock wave begins its continuous expansion, implying that there is 
still no fluid element with a positive total energy and radial velocity at that time, and that 
it is only after another ~ 100 ms that some of the fluid elements gain enough energy to be 
ejected. 

4. Conclusion 

We have done two-dimensional simulations of SASI in the supernova core to study its 
linear growth and nonlinear saturation in detail together with its role in the shock revival. 
We employed the realistic EOS and neutrino interaction rates for emission and absorption 
on nucleons. The systematic numerical experiments have been conducted, varying neutrino 
luminosities from the central object and solving spherically symmetric steady accretion flows 
through the stalled shock wave as unperturbed initial conditions. We have studied the models 
both with and without a negative entropy-gradient in the initial configurations. We have 
done the mode analysis based on the spherical harmonics decomposition. Not only the i = 1 
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single-mode velocity perturbation with an amplitude of 1% but also the random multi-mode 
velocity perturbation with an amplitude < 1% everywhere has been investigated. 

We have demonstrated that the non-radial instability grows, indeed, exponentially in 
all the models. In a model with a relatively high neutrino luminosity, we have found a 
continuous expansion of the shock wave, which will eventually lead to explosion. Since this 
model did not produce explosion in the spherically symmetric simulations, we can infer that 
the instability tends to lower the critical neutrino luminosity for the shock revival. 

In the mode analysis for the models with the initial i — 1 single-mode velocity perturba- 
tion, we have clearly seen both the linear and nonlinear growths of the instability. We have 
obtained the growth rates and characteristic oscillation frequencies by fitting the numerical 
data and found that the oscillation frequency reflects the period of the cycle of the inward 
advections of the velocity- and entropy-fluctuations and the outward propagations of pres- 
sure fluctuations. This fact suggests that the instability is indeed driven by such cycles. This 
is, however, in sharp contrast with the recent findings by Blondin & Mezzacappa (2005) that 
the repeated propagations of pressure fiuctuations are solely responsible for the instability. 
This apparent discrepancy may be ascribed to the difference of the treatments of neutrino 
coohng. 

We have observed for the models with the random multi-mode velocity fiuctuation that 
the linear phase lasts for ~ 100 ms, in which the modes with i — 1,2 are dominant. It is also 
interesting that in the nonlinear phase, where the nonlinear couplings generate all modes 
and the saturation occurs, the modes with £ — 1,2 are still dominant. Detailed analyses of 
the mode couplings will be necessary to reveal the physics behind this result. In the models 
with a negative entropy-gradient, the broadening of the spectra has been observed, which 
we have ascribed to the convective motions. This appears to be supported by the recent 
analysis of the convection in the accretion flows in a supernova core by Foglizzo et al. (2005). 
Applying their results to our models, we have shown that the convection plays a minor role 
in our models. The relative significance of SASI and convection, however, depends on the 
situations and SASI appears to be dominant for relatively low neutrino luminosities (Scheck 
et al. 2004). It is important to note here that the growth of SASI is rather slow and takes 
~ 100 ms. Although we have found a shock revival due to the SASI-enhanced neutrino 
heatings in some cases as mentioned above, whether we find such a revival in the reality 
crucially depends on if the neutrino luminosity can be sustained for this period, not an easy 
job. 

If SASI alone cannot get the job done, we will have to find something more to boost the 
shock revival. We are currently working on some combinations of microphysics and SASI 
and will soon publish the results elsewhere (Ohnishi et al. 2005). The influence of rotation 
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on SASI should be addressed also and may be important for the possible correlation of the 
kick velocity and the rotation axis (Wang et al. 2005; Ohnishi et al. 2005). 
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L^ = 6.0x10=^ erg/s 
Ly = 5.5x1 0„ erg/s 
Ly = 3.0x10^ erg/s 
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Fig. 1. — Physical quantities in unperturcbcd spherically symmetric accretion flows. Solid, 
dashed, and dotted lines represent the solutions for Lj, = 6.0 ■ 10^^, 5.5 ■ 10^^, and 3.0 ■ 10^^ 
erg/s, respectively. 
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Fig. 2. — Temporal evolutions of shock radius in the spherically symmetric simulations for 
the models given in Fig. 1 together with for L^, — 6.5, 7.0 • 10^^ erg/s. The relative deviation 
from the initial value is plotted 
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Fig. 3. — Entropy- (the left half of each panel) and density- (the right half) distributions in 
the meridian section for 1% of the i = 1 single-mode velocity perturbation. L^, = 5.5 ■ 10^^ 
erg/s is assumed for the left panels and L^, = 6.0 ■ 10^^ erg/s is for the right panels. 
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Fig. 4. — Temporal evolutions of the i = component in the spherical harmonic decompo- 
sitions. The relative deviation from the initial value is plotted. 
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Fig. 5. — Temporal evolutions of the normalized amplitudes of the £ — 1,2 modes for the 
model with L,, — 5.5 • 10^^ erg/s. The dot-dashed line represents the fitting in the hnear 
phase. 
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Fig. 6. — Temporal evolutions of the spectra in the spherical harmonics decomposition for 
the models with L^j = 3.0 ■ 10^^ erg/s (left panel) and with Lj, = 5.5 ■ 10^^ erg/s (right panel). 
The random multi-mode velocity perturbations are initially added. 
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Fig. 7. — Net heating rates in the meridian section for the models with L^, = 5.5 ■ 10^^ erg/s 
(left panel) and Lj, = 6.0 ■ 10^^ erg/s (right panel). 1% of the i = 1 single-mode velocity 
perturbation is initially added. The left (right) half of each panel represents the result for 
100 (200) ms. 
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Fig. 8. — Temporal evolutions of the net heating rate integrated over the region inside the 
shock wave. The solid and dashed lines represent the models with = 6.0 ■ 10^^ erg/s and 
= 5.5 • 10^^ erg/s, respectively. 1% of the £ — 1 single-mode velocity perturbation is 
initially added. 
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Fig. 9. — Temporal evolutions of various energies integrated over the region inside the shock 
wave for the model with L^, = 6.0 ■ 10^^ erg/s. 1% of the i = 1 single-mode velocity 
perturbation is initially added. See the text for the definition of the explosion energy. 
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Table 1. Key Variables in SASI. 



(1052 erg/s) 


7 (s-l) 


u; (s-i) 


■Rs.equil (10^ cm) 


Ws (10^ cm) 


t^adv (s ^) 




t^cyc (s ^) 


3.0 


42.4 


915 


66.9 


30.9 


1098 


5867 


925 


5.5 


45.7 


277 


128 


79.2 


262 


1832 


229 


6.0 


38.3 


188 


144 


93.5 


207 


1497 


181 


6.5 


35.6 


143 


167 


114 


159 


1199 


141 



Note. — Li, represents the model luminosities. The growth rate and the oscillation frequency denoted as 7 and w, 
respectively, are obtained by the least square fitting to the numerical results in the linear regime. Rg.equii is the initial 
shock radius and Ws is the distance between the shock radius and the neutrino sphere; Wa = Ra.oquil — The frequencies 
associated with the advection and the sound-propagation between the shock and the neutrino sphere are denoted as Wadv 
and Wsndi respectively, and are defined as Wadv = 27r/ /^^(l/iir)(if' and Wgnd = 27r//^^'' (l/cs)cir, respectively. They are 
evaluated numerically for the initial conditions. The characteristic frequency of SASI is given by the cycle frequency, 
oJcyc = 27r/[/^°(l/iir)(ir + J^" {l/cs)dr]. See the text for more detail. 



